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NUMERICAL MODELING OP POLLUTANT TRANSPORT USING A 


LAGRANGIAN MARKER PARTICLE TECHNIQUE 

By Malcolm Spaulding 
University of Rhode Island.* 

SUMMARY 

A dorivation and code have been developed for the three-dimensional mass 
transpo. V w^uation, using a particle-in-cell solution technique, to solve 
coastal zone waste discharge problems where particles are a major component of 
the waste. Improvements in the particle movement techniques have been suggested 
and typical examples illustrated. Preliminary model comparisons with analytic 
solutions for an instantaneous point release in a uniform flow have shown good 
results in resolving the waste motion. The findings to date indicate that this 
computational model will provide a useful technique to study the motion of sedi- 
ment, dredged spoils, and other particulate waste loads commonly deposited in 
coastal waters. 

INTRODUCTION 

The finite difference and finite element techniques have enjoyed widespread 
use in solving the advection-diffusion equation for water pollution transport 
problems. This fact is well documented in the bibliography on numerical models 
by Spaulding and Gordon In the application of these techniques to ocean 

waste disposal (dredged spoils or sewage sludge) and sediment transport 
several marked shortcomings are noted. Tl\e most serious problems 
encountered are: (1) generation of fictitious diffusion; (2) severe depression 

of concentration fields upstream of point discharges; (3) poor representation 
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of constituent movoinont when waste consists cl varying particle sizes (each 
with its own settling characteristics); (4) inability to readily handle 
resusponsjon of waste material on a mass dependent basis;, and (5) extreme 
difficulty in handling scale dependent diffusion processes. 

In an effort to develop a new technique which would account for these 
problem areas, the literature on solutions to the advective-diffusion equation 
was reviewed. Following recommendations made by Johnson in a review of 
ocean dredge waste disposal modeling for the Corp. of Engineers and Lange et al. 
C3,4) their work on air pollution problems, a particle-in-cell method appears 
particularly promising. This method has received considerable development in 
application to air pollution problems as noted in the work of Sklarew et al. 

, Hirt and Cook and Lange et al, The investigations have shown 

the method capable of performing useful computations to determine the fate of 
pollutant constituents under very complex air flow conditions. 

Employing the experience gained in the air pollution field, it is felt 
that the particle-in-cell method provides a very useful algorithm for handling 
complex particulate water pollutant problems. 

SYMBOLS 

C concentration of a dissolved or suspended constituent C>ng/1) 

t time 

advective velocity vector + V^J + 

^ del or gradient operator ^ ^ ^'^9y ^ ^ ^ ^ 

diffusivity tensor 

Ujj diffusion velocity vector (Dpi + V^T + W^K) 

U total pseudo transport velocity U = U . + ij 

1 1 A U 


2 


Ax 

Ay 

Az 

x,y,z 

At 

^new 


old 


T particle 


XP 

YP 

V 

M. 

1 

L 

A, 


FD 


^exact 

a 

Q 

U,V,W 


a ,a ,0 
x' y X 

D 

mean 


finite difference representation of C located at iAxjAy,kAz 
finite difference grid spacing in the x direction 

finite difference grid spacing in the y direction 

finite difference grid spacing in the z direction 

f 

three dimensional Cartesian coordinate axes 
time increment 

diagonal or principal terms of the diffusivity tensor K. . in the 

1 j 

x,y, and z directions, respectively 

updated position vector of a particle = xi yj + zk) at 

time t + At 

previous position vector of a particle at time t 
total transport velocity interpolated to the particular particle 
location 

X position of a particle within a cell referenced to the cell wall 
y position of a particle within a cell referenced to the cell wall 
cell or grid volume (AX) (AY) (AZ) 
particle mass 

characteristic dye patch size 
dissipation parameter 

finite difference approximation to the diffusion velocity 

exact differential expression for the diffusivity velocity 

2 2 2 

standard deviation o =a +o“+o 

X y y 

quantity of waste released 

x,y, and z directed velocities, respectively 

standard deviations in the x,y, and z directions, respectively 

deviation of mean from the desired mean 



a number of standard deviations for a specified confidence interval 

number of particles 

*^xo'®yo'^zo standard deviation in the x,y, and z directions, 
respectively, at t = 0. 

COMPUTATIONAL MODEL - WAPIC 


WAPIC (water advection particle-in-ccll) is a three-dimensional numerical 
solution code to the time dependent advection-diffusion equation for water 
pollutant transport employing the particle-in-cell technique. Pollutant 
concentration is represented statistically by embedded Lagrangian marker parti- 
cles in an Eulerian grid. The discussion which follows will outline the 
details of the method to include its procedural steps, stability and accuracy, 
and program options. 

Pseudo Velocity Approach 

Thti pseudo velocity approach consists of the following. Given the transport 
diffusion equation 


*K. .^C '■ 

3t A ij 

where C is the concentration of the waste, K. . tlie diffusivitv tensor, and 

X3 

lI^ the advection velocity (from a suitable hydrodynamic model or mass 
consistent velocity field). Rewriting Eq. (1) in a flux conservative form with 
the incompressibility assumption yields; 


9t 


ecu, - 


K. . 


C 

If one defines a diffusion velocity as 

^ ^ 




= 0 


( 2 ) 

( 3 ) 


4 



then a total velocity, including b iG advection velocity and the diffusion 
velocity can be written as 







Using Eqs, (2), (3), and (4) the resulting advoctivc-diffusion equation , jmes . 


at 


+ 



0 


(5) 


Equation (5) forms the fundamental starting point on which the method is 
based. Tlio solution to this equation is performed in two steps: an Eulcrian 
step and a Lagrangian step. Each will be detailed in the following paragraphs: 
Eulerian Step .- From a space staggered grid system as shown in figure 1, 
the concentrations, C, are used to calculate the diffusivity velocities 
(relationship defined in Eq, (3)). A typical finite-difference representation 
might be 


2K 

D (C. . . + C. “ , u) 
i,j ,k 1 + l,j .k-* 


C. , . C. . . 
1 + 1.3 .k 1,3 .k 

Ax 


(63 


for the X velocity component. It has been assumed that K. . reduces to K. . 

1 J IX 

such that = K . In general, it is not necessary to make this approximation, 
XIX 

however, it is commonly done. The diffusion velocities for the other coordinate 
directions are similarly defined as: 


'"i.3 + 1 ,k ~ ‘"i,! ,k 

D ■ " (C, . + C. . ^ ) Ay 
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ZKz 


D (C. 


i,j ,k + 1 i»3 ,k 


^i j k + 1 ^i.i ,k 

+ C, , . 3 Az 
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with these difference approximations the diffusion velocities are seen to be 
defined at the same points in the grid system as the advective velocities. This 
scheme proves to be very useful when calculating the total velocity at any 
position in the grid system. c 



The final task to be porforinotl in this step of the computation is the 
addition of the advective and diffusion velocities to determine the total 
or pseudo transport velocity. Since both velocities have been defined at the 
same location in the grid system, a simple addition for each coordinate direction 
velocity is all that is required. 

Lagrangian Step .- Each marker particle contained in an Eulerian cell is 
transported for one time step At accoi'ding to the total transport velocity 
U,p (defined in Eq. (4)) 

^new “ \ld ^ S’ particle (9) 

where is the Lagrangian particle position at time step t, is the 

position at time step t + At and is thr bilinear weighted 

velocity at the particle for time step t. Figure 2 shows a simple two- 
dimensional example for detouiiining the velocities at a particular particle 
location. Tliis approach can readily be generalized to the three-dimensional case. 

Although the bilinear weighting scheme or a light variation has been used 
to determine the velocity at particle locations in the models of Lange and 
Knox for the air pollutant transport problems, it is felt that this tech- 
nique is not adequate. In particular, when the streamlines or pathlines of the 
flow field are curved, the marker particles tend to drift in such a manner that 
the true streamlines and those given by the marker particles are not coincident, 

Sy reducing the time step At the error can be reduced substantially; however, 

.It the expense of computational time and, therefore, cost, 

(73 

Forester has proposed an iterative approach to particle motion 
•orediction that alleviates many of the problems inherent with the simple 
bilinear velocity weighting scheme. Instead of employing Eq, (9) to represent 
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the particle motion computation, he suggests the following; 

I) - 4 - ( it + It \ 

Slew old It particle T particle} *T 


( 10 ) 


whore IJ,^, particle the bilinear weighted velocity at the old 


»- r new 

particle position and U,^ represents the bilinear weighted 

velocity at the new particle position With an iteration procedure where 


->• 


each estimate of i® calculate an updated particle velocity 

particle solved to find the new particle positions , for a 

given convergence error on 

In an attempt to assess the usefulness of the iteration procedure, several 
trial cases have been executed, Using a simple four-cell velocity field with 
rotational currents (ns shown in figures 3 and 4) , a grid spacing of 10 m in 
both tlie X and y directions, and varying time step, particle trajectories 
were predicted with and without the iteration procedure previously outlined. 
Figure 3 shows that without iteration and even for very small time steps Ci,e., 
1/10 of a grid space per time step) the trajectories show a definite spiraling 
outward from the true particle pathline--a closed circle for this velocity field, 
Figure 4 presents the same simulations but employing the iteration procedure. 


It is readily seen that there is a marked improvement in particle trajectory 
representation with this approach. 

Therefore, it is concluded that the use of the proposed iteration tech- 
nique is the best method to determine particle motion, since It requires less 
computational effort for more accurate results. Without this procedure, it can 
be seen that eddy type structures are more rapidly depleted of their pollutant 
leads than occurs in the real case. 

The final portion of the Lagrangian step is to compute the concentration 
given the new particle positions and mass. A "fictitious" volume is associated 


7 



with ouch ptirticio mass at the start of the computation. The ovorlap of this 
volume into the surrounding hulerif-n colls do'cx*mincs tlio distribution of mass 
into the various colls and allows the concentration field to bo calculated. 
Figure 5 shows the details of the computational procedure for a typical particle 
in a two-dimensional grid system. This algorithm is simply repeated for all 
particles in the system for a given time. A further refinement can be added in 
the Lagrangian step to allow for gx’id expansion or a grid system that moves 
with the mean motion of the particles. Incox’poration of this feature permits 
maximum resolution of the concentration field with a minimum number of computa- 
tional grids, and is particularly useful for determining the behavior of "puff" 
or instantaneous point releases. Variations of the moving grid system procedure 
would allow tho resolution of some particular portion of the pollutant field, 
i.e,, determining where particles of a spocified mass are transported. 

Figure 6 presents a flow chart for the computational procedure that has 
been outlined in the previous paragraphs. 

Velocity Information 

In order to perforin computations with IVAPIC, a mass consistent velocity 
field is required. That is a velocity field that conserves water mass. This 
information for both spatial and temporal variations can be obtained from a 
finite difference or finite element model of the coastal cone hydx'odynamics . 
Reference 1 gives an extensive list of suitable circulation models of both the 
two-dimensional vertically averaged and three-dimensional types. It is also 
possible to use statistical predictor type models based on data taken for 

a particular area Tiie main requirement of each, however, is conservation of 
water mass. 
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Tlie final rcciuiroment of the flow field is that the velocity Information 
can be resolved into the space staggered grid system ns shown in figure 1. A 
simple linear weight processing technique witli a water conservation constraint 
should readily achieve this goal. 

Boundary Conditions 

WAPIC has two basic typos of boundary conditions, a closed or zero mass flux 
boundary (U.j. pgrticlos*^^ = 0, and a mass flux boundary pa^tlcles^^ “ constant, 
In the application of these boundary specifications, consistency of boundary 
type for both the Eulerian and LagraAgian step calculations as well as the input 
velocity field must be maintained. Specifically, the Eulerian step must account 
for the flux of pollutant particles due to diffusion while the Lagrangian step 
resolves the advective particle flux. 

Another boundary condition particularly applicable to the sediment trans- 
port problem is a deposition or storage boundary. When particles which are 
settling from the flow field reach the bottom of the water column they may bo 
stored at that geographic location until such time as they might be resuspended 
by an increased shear at the water-soil interface. 

Diffusion Parameters 

In the general case, WAPIC can, in principle, accommodate the full Cartesian 
eddy diffusivity tensor K. . . In practice, isotropy of the water turbulence is 
assumed and only the diagonal terms noted as K , K . and K_ remain. In 

X > - 

addition, for most coastal zone flow conditions it is generally assumed 
It still remains, however, to determine values or empirical relations that 
approximate the diffusion processes for coastal cone areas. 
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Tlie majority of iiivostlgntions for horizontal turbulent diffusion in tho 
ocoun monitor tho size of tracer dye patch as a function of time and then using 
si~vU diffusion models, which disregard both shear current and vortical trans- 
port, calculate tho diffusion coefficients 

The values of horizontal diffusion coefficients obtained In tlio ocean range 
2 8 2 

from 5 X 10 to 4 X 10 cm /soc. Tho largest number of measurements were 
obtained at tho ocean surface, In gonoral, it was found that the was a 
direct function of the diffusing patch size L-, Figure 7 shows a collection of 
data and illustrates this relationship quite well. Although tho data scatter 
is significant, in general 

s! A, 10 “^ft < b 5l0^ft 

X u 

where Aj^ is a constant called the dissipation parameter and varies from 0,005 
to 0.00015 ft^'^^/soc. 

Diffusion coefficients obtained in the ocean are not strictly applicable 
to estuarine waters? they are, however, useful in demonstrating trends and 
magnitudes. Brandsima and Divoky have shown that data from tidal estuaries 
lie right on tho ocean data when applying tho four-thirds power law. Tliere- 
fore, Eq. (11) should be a useful first approximation to horizontal diffusion 
parameters for both coastal zone and estuarine waters. Further refinements to 
hori 2 ont.ii diffusion prediction will probably rely on field experiments in the 
local area. 

Tiie vertical diffusion coefficient in the ocean is generally much smaller 
than the horizontal coefficient, Table I shows a summary of typical values 
employed for a variety of coastal and open ocean conditions. Although there 
appear to be no general relations for vertical diffusion, it normally displays 
a maximum near the surface (caused by wind mixing) and decreases with depth. 
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In addition, the vertical mixing decreases us the Richardson number (measure of 
density stratification) increases. Table IX Indicates some typical formulations 
for the vertical diffusion coefficient. 


Accuracy and Stability Requirements 

In all numerical solutions to tlie advectivo-diffusion transport equations, 
there are requirements which must bo met in order to achieve accurate and stable 
solutions to the equations. Discussed below are detailed conditions for the 
pseudo velocity approach. 

Number of Eulorian Cells per Problem ," It is a basic requirement for all 
methods that use Eulerian grid systems and finite difference or elements 
opproximatioiif sufficient grids are used to resolve the concentration 

gradients of interest. When this condition is achieved, the numerically 
calculated gradients of concentration and t)ie actual gradients agree with one 
another. In addition, as the cells used to represent the concentration field 
increase the agreement between actual concentrations and those of the numerical 
computation becomes better. 

In an effort to develop an approximate quantitative estimate of the lower 
limit and the number of grids necessary to resolve the concentration distribu- 
tion, we will follow the approach outlined by Lange . Consider the one- 

dimensional finite difference approximation to the diffusion velocity given by.‘ 


K 


+ 1/2 


^i + 1/2 + 1 " 

Ax 


U21 


C. 


"i + 1/2 

Using a Taylor series expansion about i + 1/2 for + 2 ^ ^ibd C^, substitut- 
ing into Eq. 12 and dropping the i + 1/2 subscript yields 


U = - 


AxC 


1 T 

Ax"^ + higher order terms 

"" i 


(13) 
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Tliou expanding Eq. (13) and dropping higher order terms gives 

K.. 

U = ~ 


\ ^ 

C ■ 24 C 




CH3 


3x‘ 


The first term on the right-hand side of Eq.(14)is the exact differential 
expression for the diffusivity velocity. Taking the ratio of the diffusion 
velocity defined by finite differences (Oq. (14)) and the exact differential 
results in: 


U 


FD 


3x^ Ax^ 


U 


exact 


3C 

3x 


24 


(IS) 


In order to obtain a quantitative measure of this ratio estimates of tlie 
first and third derivatives of concentration with distance are required. 
Assuming a simple one-dimensional Gaussian die. fribution: 


= a 


xV2a^ 


C = ^ e 
a 


and finding the first 


(16) 


and third deri’fatives 




x^/2a^ 


3x 


3JC_ „ gx 


2 i 2 2 

, xV - x /2a^ 

3 - ® 

. o j: 


ax'" a- 

this estimate can be made. 

Substituting Eqs,(17) and (18) into Eq. (15) one obtains 


(17) 


(18) 


U 


FD 


exact 


AC 


3| 4 
7 


(19) 


Employing the fact that 99.9 percent of all particles in a one-dimensional 
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(Gaussian distribution are within x “ 3<J, we can substitute x « 3or in 


;Gq, (19) and find 



^exact 


1 



( 20 ) 


Prom this relationship, it is seen that when a grid cell is either larger or 
comparable to the standard deviation for particle distribution, the finite 
difference algorithm underestimates the diffusion velocity. However, as more 
grid cells are employed for a given standard deviation, the finite difference 
approximation becomes an increasingly better estimate of the diffusion velocity, 
.Number of Particles per Cell ,- Since the particles represent a statistically 
quantized density, it is desirable to have as many particles as possible for any 
particular problem. An upper bound on ppr^’cles would be either the storage 
capability or maximum program run time (this is directly proportional to the 
number of particles) for a particular computer. The lower bound is at least 
one particle per cell. Fewer particles and in particular when a particle has 
no neighbors withi/i one cell length, the particle is moved to the grid boundary 
and "frozen" in place when considering the diffusion velocity. This freezing 
process is a direct consequence of the algorithm for determining concentrations 
from particle positions (see figure 5). It can be readily seen that any parti- 
cle must distribute its mass to its nearest neighbor cells. If another distri- 


bution algorithm were chosen which distributed the particles over a greater 
volume, the one particle per cell minimum could be relaxed. 

Particle Generation .- In general, it is desirable to start either steady 
state or time dependent point discharge simulations with some initial Gaussian 
or normal distribution. Several experiments were conducted using random number 
generator techniques similar to those in the work of Lange , which proved 
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to give distributions that displayed noticeable '’bunching" of particles around 
the directions of the coordinate axes. It was felt that this problem was 
directly related to the quality of the normal random number generator available 
at the time. 

In an effort to alleviate this problem, it was decided to specify the exact 
number of particles in each cell using the Gaussian probability distribution 
function. Knowing the cell size, standard deviation of the distribution, and 
the total number of particles in the field one can then pompute the number of 
particles in each cell using the normal distribution constraint. Next, a 
uniform random number generator was used to locate the appropriate number of 
particles in each cell. All cells were filled in this manner. If it is 
decided to obtain another type of initial condition, all one needs to do is 
specify the desired distribution as a constraint for defining the manner in 
which these cells are to be filled with particles. The remainder of the algorit.im 
remains unchanged, 

Time Step Restriction .- Similar to many finite difference schemes, hAPIC 
has a restriction that the fastest moving particle cannot move more than one- 
half cell length in one time stop. At. In equation form this restriction 
becomes 

fU At V At WJ^t ■) . 

maximum > _! > { 1 

L Ax Ay Az J (21) 

where U.j,V,p and W„ are the total or pseudeo transport velocities At is the 
time step, and Ax, Ay, and Az are the grid spacings in the x,y, and z 
directions, respectively. These are accuracy conditions rather than stability 
ones. 

As was discussed earlier, if one does not choose to employ the iteration 
particle location technique and the flow-field displays marked eddy structure 
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or curved streamlines, it becomes necessary to further restrict the maximum time 
step. The exact value to be employed in that particular situation is highly 
problem dependent and it is, therefore, recommended to use the iteration 
procedure in all cases. 

.Number of Particles for Given Confidence Limit on Mean. - When undertaking 
a problem, it is extremely desirable to have at least a first estimate of the 
number of particles necessary to define a particular distribution or starting 
condition. Since the Gaussian distribution is often employed to define the 
initial conditions, a logical manner to determine the number of particles is to 
use simple statistical the''ry to find the number of particles necessary to 
assume a given confidence in the deviation of the mean from the desired mean. 

In Bowker and Lieberman's text it is shown that the lower one-sided 

deviation from the mean confidence limit is; 


D = K 
mean a 




1/2 




/iT 


2 2 ’2 2 
a = + O +0 

X y z 

where D - deviation of the mean from the desired value, a - the standard 
mean 

deviation of the normal distribution, n-the number of particles, and — the 

number of standard deviations for a given confidence Interval (e.g., for 95 

percent confidence interval - K„ = 1.645). 

In an effort to gain a more quantitative estimate of this relation, a 

parameteric study on confidence limit (70 - 95 percent) with a a = a 1000. m, 

X y 

= 1.5 m, and a mean of 5000 m is shown in figure 8. In general, it indicates 
that for this distribution about lOOO particles are adequate to define the 
initial conditions . 

The influence of a on the number of particles for a given confidence 
limit and deviation from the mean is shown in figure 9. As expected the 
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smaller the standard deviation the loss particles necessary to define the mean 
to within a given confidence i> erval , 

Although figures 8 and 9 are for a particular problem* the technique out- 
lined provides a simple procedure to determine an order of magnitude estimate 
of the number of particles necessary for a given problem. In addition, pre- 
liminary estimates on error bounds for subsequent program solution can be 
obtained. 

Options in WAPld 

Sources and Sinks .- Sources or sinks fy. pollutant can be generated anywher 
within the Eulerian grid mesh and may eit ter be instantaneous, continuous, or 
intermittent. From accuracy considerations, previously discussed, the 
initial pollutant distribution must cover at least two grid lengths, otherwise, 
the particles diffuse too slowly. This, therefore, requires special sub-grid 
treatment if accurate near-field results are necessary. 

Deposition and Resuspension . - The deposition of particles can readily be 
laccommodated with IVAPIC, Knowing the particle density and shape/surface 
characteristics, a settling velocity can be approximated for each particle and 
simply added to the total transport velocity, The ability to handle each 
particle separately allows highly detailed behavior of waste material contain- 
ing varying size distribution such as sediment or dredged material to be 
predicted, 

Resuspension of material from the bottom or some other storage area can be 
handled on a particle size basis provided the hydrodynamic model input provides 
sufficient detail of bottom current structure. 

Flocculation . - A rather prevalent process in the settling of fine sediment 
and waste material is flocculation. Electrostatic forces between neigliboring 
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particles cause an aggrograto to be formed with new settling clmracteristics . 
Classical finite difference and finite element solutions to the advoctive- 
diffusion equation arc unable to represent this process very accurately since 
they rely on a gross adjustment, such as an additional sink of the waste 
material to replicate the flocculation process. 

WAPIC allows one to assign electrostatic and other properties to a given 
particle and the separation distance between particles is easily obtained 
from the Lagrangian marker paths. With this information, better resolution 
of the details of flocculation can be achieved. 

Time History of Particles. - In some sediment transport studies, radioactive 
tracer materials are used to monitor the time dependent sediment motion. Using 
the detailed time history of the particles available with WAPIC, one can 
assume given radioactive decay rates for the tracer particles as they travel 
along their trajectories. These time path histories, therefore, provide an 
additional capability over classical solutions to the advective-diffusion 
equation . 

VERIFICATION 

In an effort to perform preliminary model verification, several simple 
instantaneous releases of pollutant were simulated. A spherically symmetric 
Gaussian puff distribution consisting of 1984 particles was instantaneously 
released into a uniform velocity field. The model parameters for this case 
were : 

Grid spacing Ax = Ay = 500 m Az = 1.5 m 

2 2 

Diffusion Coefficients K = K ,= 10 m /sec,, K = 0.0001 m /sec. 

^ y 2 
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Standard Deviations 

(of initial distributions) = 1000 ° 

Time step At = 250 see. 

Grid system with mean particle motion, 

A corresponding analytical solution for this problem was obtained from 
the work of Ikubo ot al . and modified to include an initial distribution 

of the concentration field. The solution is given by 


C(x,y, z,t) = 


(2tt) 
X exp 


3/2 i 


xo ♦ 


Ko' * 


1/2 


I Cx - ut)^ 

^ a + 2K t 
xO X 


a + 2K 
yO y 


—4 — )i 

< t 0 .“ + 2K t / 
y zO z • j 


(23) 


Comparison of the model results with the analytic solution are sho^^m in 
figures 10 and 11 . Figure 10 gives a comparison of the concentration distri- 
bution in downstream direction for several time increments. Only half of the 
profiles were plotted on the first and last time steps in order to aid in 
interpreting the graph, 'Hie variation of model predictions to data is with 
+ 5 percent for all cases, Figure 11 shows a similar graph but in the vertical 
plane. Again, the comparison appears within 5 percent. If an increased number 
of particles were employed and better initial resolution (finer than 


~ °^y0 “ ^zO " error in model predictions could be further reduced. 

CONCbUSIONS 


A derivation and code have been developed for the three-dimensional mass 
transport equation, using a particle-in-cell solution technique, to solve 
coastal zone waste discharge problems where particles are a major component 
of the waste. Improvements in the particle movement techniques have been 
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suggnsted and typical examples illustrated. Preliminary model comparisons with 
wiaJytic solutions for an instantaneous point release in a uniform flow have 
shown good results in resolving the waste motion. The findings to date indicate 
that this computational model will provide a useful technique to study the 
motion of sediment, dredged spoils, and other particulate waste loads commonly 
deposited in coastal waters. 
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TABLE I 


SUMMARY OF VALUES 

OF VERTICAL DIFFUSION COEFFICIENT K» IN THE OCEAN 
(Reference 10) ^ 

Note: Molecular diffusivity for heat: I.S x 10* ** ^ cm^/iec (at ZO^^C, I atm) 

■alt: 1.3 X 10*^ cm^/aec (at 20*^C, 1 atm) 




Vertical Diffusion 

1 

1 

Current or oceanic 

Depth of 

Coefficient K, 
(cm^/*ec) ^ 


region 

layer (m) 

Reference 


Philippine Trench 

5000-9788 

Z.0-3. 2 

Schmidt, 1917 

Algerian Coast 

0- 20 

35-40 

Schmidt, 1917 

Mediterranean 

0- 28 

42 

Schmidt, 1917 

California Current 

0- 200 

30-40 

McEwen, 1919 

Caspian Sea 

0- 100 

1-3 

Stockman, 1936 

Barents Sea 


4-14 

Subov, 1938 

Bay of Biscay 
Equatorial Atlantic 

O 

o 

1 

o 

Z* 

Fjeldstad, 1933 

Ocean 

0- 50 

320 

Defant, 1932 

Randesfjord 

0- 15 

0. 1-0.4 

Jacobsen, 1913 

Schultz Grand 

0- 25 

0, 04-0. 74 

Jacobsen, 1913 

Kuroshlo 

0- 200 

30-80 

Sverdrup-Staff, 

Kuroshlo 
Southern Atlantic 

0- 400 

7-90 

Suda, 1936 

Ocean 

400-1400 

5-10 

Defant, 1936 

Arctic Ocean 

200- 500 

20-5C 

Sverdrup, 1933 

Carribean Sea 

500- 700 

2.8 

Seiwell, 1938 

South Atlantic Ocean 

3000-Bottom 

4 

Defant, 1936 

South Atlantic Ocean 

Near Bottom 

4 

Wattenbere, 19] 


West Atlantic Trough 

(50°Sto 10° N) 

North Atlantic 
Indian Ocean 

Pacific Ocean 

Tidal Channel 

(Mersey estuary 

and Irish Sea) 

Near Cape Kennedy, 
Florida 
Bikini Lagoon 

Coast of benmark 
California Coast 


Near Bottom 


Near Bottom 

0- ZO 
(bottom) 

Surface Layer 
0- 50 

(bottom) 


2-40 

19 (in August) 
1. 3 (in Summer 
260 

0.05-1 
0. l-IO 

15-180 

(at wind force 
3-4) 


* As given by Defant, 1961 

** As given by Bowden, 1962 
As given by Harrcmocs, 1967 
As given by Wiegel. 1964 


WSst. 1955 


Koezv, 1956 *■<' 

Bowden, 1965 
(with R, from 
0. 1 to ^ 2. 01 
Carter and Okubo, 

1965 

Monk, Ewing and 

Rcvelie, 1°49 

Harremoesi )^o 7 

Foxworthy, Tib / 

and Barsom, 1°66 

St.’mmel and 
Woodcock, 1951 
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TABLi: II 


SUMMARY OF FORMULAS 

ON 'CORRELATION OF VERTICAL DIFFUSION COEFFICIENT 

WITH RICHARDSON'S NUMBER R, (OR DENSITY GRADIENT 0) 
(Roforonco 10) ^ 

Note: K Q : K at R. = 0, l.e,, the neutral case fl : proportionality constant varlea 
* ^ from case to case 


Rossby and Montgomery 


& 


(1935)X' 



Rossby and Montgomery 


S 

+B «,)•” 

(I935)x« 



Holzman {1943)* 

K 

z 

tf- 


Yamamoto (1959)’*' 

K 

z 

r 


Mamayev (1958)’*' 


= 


Mtink and Anderson 

K 



(194 8)X'X' 

z 





=3.33 based upon data by Jacobsen (1913) 




and Taylor (1931) 

Harremocs (I960) 


s 

- ,rt-3 -2/3 2/ 

5x10 X c cm /sec 

note! G lnm“S approximate experimental 

- 9 - 5 - 1 

range 5x10 <c< 15x10 m 

Kolesnikov, et al 
(1961)x»ti>.v 

Kz 


K, , + — In cm^/sec 

^ mln c 

tnin ® empirically determined 

to be: 

K- , =12, e = 8.3 X 10'^ (1958 and 

1960 observations) 

K , = 2, e = 10.0 X 10“^ (1959 

2 observations) 

Koh and Fan (1969) 


- 

10 /C (K^ in cm^/sec; c In m”^^ 




4 X lO""^^ e lO-^m"^ 


* As given by Okubo 
As given by nowclcn (I 0 (j 2 ) 

The formulas presented in the translated version arc apparently erroneous, 












Position Coordinate (M) 






Position Coordinate CM) 


/ y / 


Velocity - 1 M/sec. 

Grid Spacing 

X AX - 10 M 

Y AY - 10 M 

Time Step AT - 1, and 5 Se 
Convergence Error - „01 


. , 5. 


/X ^ ' / y /■ y y y 

15 20 25 

X Position Coordinate 


Grid System Showing Velocity 
Field 


Figure 4. - Particle Trajectories for Varying Time Steps (1 and 5 Sec,) with 
Position Iteration, 






CONCENTItATION COMPUTATIONS 

^I,J = ^ CiZ) - XP) - YP) 

^I+1,J = ^ im (XP - (-^ - YP) 

'^I*l,Jtl = 72 CAZKXP - % (YP - 
“2 (AZ)C^- XP}(YP - 


Figure S.- Computational Algorithm to Determine 

Concentrations, given Particle Position 
and Mass. 
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start 




Initialize Variables 




Place Particles in Cells 
to Represent Appropriate 
Loading Conditions 
(Mass and Postition) 




\ 


Read Velocity Data 
for TVrea of Interest 
(see Velocity Data Processor) 


t 


Determine Diffusion 
Velocities from Concentration 
Prediction 


Compute Settling 
Velocities Using 
Particle Mass and 
Density Difference 




Conbine Advection, Diffusion and 
Settling Velocity to Obtain Total 
Particle Velocity 


Set Boundary 
Conditions 






Move Particles in Response 
to Total Velocity Field 




Using New Particle Positions 
Calculate New Concentrations 
Dnploying Volume weighted 
Technique 


I 


C3iec)c Conservation of 
Constituent Mass 




Point/Plot Velocity 
and Concentration 
Predictions 


FIGURE 6 Generalized Flow CJ^art for WAPIC 


Increjnont 

Time 


‘-T-' 

Determine Con- 
centration Field 

in New WAPIC Grid 
System 

( 

Hove WAPIC Grid System 
in Relationship to 
Velocity Field Grid 
System Options 

a. No motion 

b. Expansion about 
a Fixed Point 

c. Move witli mean 
particle motion 

d. Combination of 
Eulcrian Lagramgian 
approach 
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Figure 7, Horizontal Diffusion Coefficient 
as a function of Horizontal Scale 
(Reference 11) 


2D 



NUMBER OF PARTICLES CN) 



PERCENT DEVIATION OF MEAN FROM DESIRED MEAN 


Figure 8*- Number o£ Particles vs, Percent Deviation of Mean for 

Varying Confidence Limits Assuming Fixed Mean (p = 5000) 
and Standard Deviation (c = 1414) . 


0 




T 


NORMAL DISTRIBUTION 

DEVIATION OF MEAN - 
CONFIDENCE LIMIT - 95% 


STANDARD DEVIATION 7 


Figure 9.- Influence of Variation in Standard Deviation of 
Number of Particles for a Given Deviation of :he 
Mean and Confidence Limit. 



CONCENTRATION 


/XNALrnC SOL'N 


o WAPIC PREDICriONS 

TAKEN AT Z=1.5 H 
ABOVE SOURCE (Z=0) 



DISTANCE DOWNSTREAM CSOOM INCREMENTS;) 

Figure 10.- Concentration \‘s. Horizontal Distance at Time Increments 0., 1000., 
2000., 3000. and 4000. For Instantaneous PT. Release in a Uniform 
Advection Field. 


CONCKNTMTION 



0 2 4 6 8 10 

VERTICAL DISTANCE (1.5 M INCREMENTS) 

Figure 11.- Concentration vs. Vertical Distance at Time 


Increments 0 and 4000 Sec. for an Instantaneous 
Point Release in a Uniform Advection Field 


